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On the relationship between instability and Lyapunov times for the 
3-body problem 
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ABSTRACT 

In this study we consider the relationship between the survival time and the Lyapunov time 
for 3-body systems. It is sho wn that the Sitnikov problem exhibits a two-part power law re- 
lationship as demonstrated in Mikkola & Tanikawa (2007|) for the general 3-body problem. 
Using an approximate Poincare map on an appropriate surface of section, we delineate escape 
regions in a domain of initial conditions and use these regions to analytically obtain a new 
functional relationship between the Lyapunov time and the survival time for the 3-body prob- 
lem. The marginal probability distributions of the Lyapunov and survival times are discussed 
and we show that the probability density function of Lyapunov times for the Sitnikov problem 
is similar to that for the general 3-body problem. 

Key words: Stellar dynamics - celestial mechanics - time. 



1 INTRODUCTION 

A correlation between the Lyapunov time, the time it takes for 
nearby orbits to diverge by e, and the time in which an orbit under- 
goes a sudden tran sition in the Solar System was first discussed by 
Lecar et al] ( 1992). From a study of orbits of asteroids I Sop er et"al] 



1990b between Jupiter and Saturn, the authors noted a relationship 
. *-h between the Lyapunov time, f/, and the time, tj, which an asteroid 
■ takes to cross the orbit of Jupiter or Saturn. A correlation between 
%— I these two ti me scales wa s also n oted for asteroids in the outer as- 
^ , teroid belt in lLecaretal.Ul992bh . In both studies it was found that 
trf and t\ are related by 



k 
c 



(1) 



where (3 is a constant, C is a normalization constant, and A is a 
constant of proportionality. 

In support of the relationship l[T), iLecar et all d 19921) consid- 
ered the elliptic restricted 3-body problem in which the massless 
particle, 1113, began its motion around the secondary mass which 
was 1/9 the mass of the primary where the orbit of the secondary 
body had an eccentricity of 0.1. In this example, was taken to be 
the time it took for /713 to escape via one of the collinear Lagrange 
points. Correlating data from 1000 orbits, the study found that QJ 
holds for (3^ 1.8. 

There have been many other investigations into the 
relationship between L yapunov times and survival times. 
iLevison & Duncan] dl993h in a study of Edgeworth-Kuiper belt 
objects showed a relationship between the Lyapunov time and 
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the time it takes for these object to cross the orbit of Neptune. 
In this study, the authors considered orbits of 200 particles with 
eccentricities between 0.01 and 0.1. They found that (QJ holds but 
with a slightly highe r exponent value (3 ~ 1.9. In another study, 
IrVlurison et al.l J 19941) considered the restricted elliptic 3-body 
problem with Jupiter as the secondary mass. Again, they found 
that the relationship |Q} holds with (3 = 1.74 ±0.03. 

Despite all the support of the relationship {JJ, there have 
been s ome disagreements with the relationship. Murray & Holmanl 
found that the relationship does not hold for some bodies in 
the outer belt. Their explanation is based on the properties of a sys- 
tem controlled by a critical KAM curve. The dynamics of the outer 
asteroid belt are not controlled by a single critical KAM curve, and 
the authors argued that this means that there is no reason to ex- 
pect a simple scaling between the Lyapunov time an d the escape 
time. In another study M orbidelli & Froeschltj |_[996 1 gave a two 
part relationship. They suggested that for orbits in the Nekhoroshev 
regime, the relationship between tj and /; should be exponential 
(i.e. t,] ~ exp(//)) whereas in a regime with resonance overlapping 
a relationship of the form iQJ can hold. 

An investigation of the relationship between escape times and 
Lyapunov ti mes for the general 3-body pr oblem has recently been 
conducted by Mikkola & T anikawal ^20071) . The authors looked for 
a correlation between the Lyapunov time and escape time for the 
planar 3-body problem. The authors considered over 10000 initial 
values in a domain of initial conditions and found that a two part 
power law works best. For orbits with small f/, the authors sug- 
gested that a power law (TJ with exponent (3 ~ 2.3 approximates 
the data whereas for large fy, the power law fits better with (3 ~ 1. 

In this study we discuss the relationship between the Lya- 
punov time and the survival time for a specific 3-body configuration 
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Figure 1. The Sitnikov Problem 



known as the Sitnikov problem. In section[2]we demonstrate that f/ 
and for orbits of the Sitnikov problem exhibit a two part power 
law relationship similar to that for the general 3-body problem. It 
is further shown that the relationship between ?/ and tj for small 
t\ is dependent on the eccentricity of the binary system. In section 
[3] we present an approximate map for the Poincare map discussed 
in iMoserl i 19731) for the Sitnikov problem. We then demonstrate 
that the relationship between t[ and t^, for orbits computed with the 
approximate map, is similar to both the Sitnikov problem and the 
general 3-body problem. Using the approximate Poincare map we 
delineate a region of initial conditions which escape quickly; orbits 
for these initial conditions are used to construct a new functional 
relationship between // and Finally, in section|4]we discuss the 
probability distributions of ?/ and and compare them to the dis- 
tributions for the general 3-body problem. 



2 SITNIKOV PROBLEM 

The Sitnikov problem is the problem of the motion of a massless 
particle, 1113, on the axis of symmetry, L, of an equal mass (mj = mq) 
binary (FigureQJ. Units are chosen such that the total mass of the 
binary is unity, the period of the binary is 271 and the gravitational 
constant G = 1 . The equation of motion for takes the form 
7 

(2) 



Vz 2 + r 2 

where z is the position of mj, along L, z = corresponds to the plane 
of the binary, and r is the distance from one of the binary particles 
to the centre of mass. The specific energy for 1113 is given by 



2 Vz 2 + r 2 



(3) 



The value of r can be computed from Kepler's equation or for small 
eccentricities, e, of the binary, we can approximate r to first order 
in e by, 



- (l-ecos(0). 



2.1 Definitions 



(4) 



The survival time of orbits for the Sitnikov Problem is defined as 
the duration of the numerical experiment to the point where 1113 es- 
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Figure 2. Survival time of orbits of the Sitnikov problem with initial con- 
ditions in the circle whose radius is given by {8}. Each initial condition is 
plotted in polar coordinates where the radial argument is determined by the 
initial velocity and the angular argument is determined by to. The colour 
associated with each point indicates the number of periods of the binary 
before escape was determined. 



capes from the system. We say m<$ has escaped at time t if sign(z(f)) 
= sign(z(f)) and 



>0. 



(5) 



where X is a con stant such that X = (1 — e)/2. It can be shown 
jUrminskvll2008bh that if the motion of m$ satisfies the above con- 
ditions, then the system's final motion is hyperbolic-elliptic. 

The Lyapunov time for orbits of the Sitnikov problem can be 
computed from the solutions of the variational equations, 8z(f). For 
chaotic systems, the magnitude of the variational solutions has or- 
der 



Mt)\ ft 

where f; is the Lyapunov time. Evaluating ((6) at time t 
solving for ?/ gives the Lyapunov time as 

t td 
' lnfl&foJI/l&ol)' 

where 8zo = §z(0), 



(6) 
trf and 

(7) 



2.2 Initial Conditions 

If we take initial conditions for {2]l such that z(frj) = for an initial 
time to = 0, we can determine from ((5} that as ms crosses the plane 
of motion of the binary, a velocity of 



(8) 



z(/oJ > \ hr 



will ensure that escapes the system without returning to the 
plane of the binary. As @ is periodic with period 271, we can con- 
sider initial conditions in polar coordinates where time is the angu- 
lar argument and z(frj) is the radial argument. Thus we can define 
the set of initial conditions for ((2) as the circle of radius \FTfk 
centred at the origin. 

Figure [2] shows the complement of the region defined by ((8}. 
A grid of initial conditions was chosen insid e the disk and iter - 
ated forward using the Bulirsch-Stoer method jPress et ail i 19921) ) 
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Figure 3. The scatter diagram of the survival time, tj, and the Lyapunov 
time, f;, for the Sitnikov problem where the eccentricity of the binary is 
£ = 0.61. The dashed line is a median curve such that at any position along 
the curve there are an equal number of scatter points above and below the 
line. 
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Figure 4. The solid curve is the median curve shown in Figure [5] The two 
dashed lines represent the power law relationship (TJ on different time in- 
tervals. For 1 < t/ < 6) we find that p as 2.5 best approximates the data. For 
// > 6 we find that p « 1 . 1 works better. 



for either 10000 time units or until the escape criterion was satis- 
fied. The colour associated with each initial condition represents 
the number of periods of the binary before either satisfied the 
escape criterion, or the numerical integration algorithm reached its 
maximum time. The outer grey region represents initial conditions 
in which the escape criterion was satisfied before the mass returned 
to the plane of motion of the binary. The inner black regions corre- 
spond to initial conditions whose orbits remain bounded. 



2.3 Results 

To demonstrate a relationship between the survival time and the 
Lyapunov time for the Sitnikov problem, 10000 initial conditions 
where chosen in the region desc ribed in section |2!2l The Bulirsch- 
Stoer method dPressetalJI 1992b was chosen as the numerical inte- 
grator with a relative tolerance of 10~ 12 . We integrated each initial 
condition simultaneously with the variational equations for 100000 
time units or until the solution satisfied the escape criterion. If the 
solution failed to satisfy the escape criterion within the time limit 
given to the integrator it was not considered in the results. This is 
because there is a large area of bounded motion, approximately the 
black regions in Figure[2] which never escape. 

Figure [3] displays th e (ti,tj) scatter diagr a m in logarithmic 
scale for e = 0.61. As in Mikkola & Tanikawa ( .2007b . the range 
of ti is divided into 50 intervals each containing an equal number 
of points. The dashed curve in Figure [3] displays the median of the 
survival ti me, tj, in each of the 50 interv als. Comparing this plot to 
Figure 3 in Mik kola & Tanik awa ( 2007) we note some similarities. 
First, for small t\ the median curve is steeper compared to larger t\ 
values. Secondly, the density of the scatter points becomes smaller 
as t\ increases. Finally, for small the scatter plot has horizontal 
band-like structures. In Figure[4]we re-plot the median curve in Fig- 
ure [3] and approximate the median curve with (JJ on two separate 
t\ intervals. For 1 < t\ < 6 we find that (JJ approximates the curve 
with p = 2.5. For t\ > 6 we find that (JJ approximates the data with 
P=l.l. 
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Figure 5. The median (7/ , tj ) curve for varying eccentricities of the Sitnikov 
problem. 



2.4 Perturbations caused by large eccentricities 

Of the systems studied by Mikkola & Tanikawa (2007), in partic- 
ular the free-fall 3-body problem, the change in energy of the es- 
caping body can vary widely depending on the interaction with the 
resulting binary system. This can be modelled in the Sitnikov prob- 
lem by increasing the eccentricity of the binary. 

Figure [5] represents the median curves for a series of exper- 
iments in which the eccentricity of the binary system in the Sit- 
nikov problem is varied. One distinguishing feature in this figure is 
the increasing prominence of a two part power law relationship be- 
tween the survival time and the Lyapunov time as the eccentricity 
increases. Interestingly, the minimum Lyapunov and survival times 
also decrease as the eccentricity increases. 

In summary, by increasing the eccentricity of the Sitnikov 
problem we can cause large perturbations to the energy of mj, and 
for e large enough the two-part power law becomes more promi- 
nent. An explanation for the different power laws between small 
and large t\ is still needed. To help provide further theoretical expla- 
nations for the relationship we can turn t o an approximate P oincare 
map for the Sitnikov problem derived in lUrminskvl ( T2008bb . 
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3 APPROXIMATE POINCARE MAP 

The plane which corresponds to z = is a natural choice for a sur- 
face of section (SOS) on which to study escape with the Sitnikov 
problem. On the SOS we can consider a map <]) : (vo,fo) ~* ( v li ; l) 
which takes m-x, from one crossing of the SOS to the next crossing. 
If m-i, is on the SOS at time to, (]) is a map wh ich brings vp = z(to) 
to time t\ > to where vi = i(t\ ) and z(t\ ) = O. lMosej 1 1973h shows 
that there exists a real analytic simple closed curve in R 2 in whose 
interior, Do, the mapping (]) is defined. In addition, § maps Do onto 
a domain D\ and for e > the boundary curves for Do and D\ 
intersect transversally. Any point not in Do is said to escape. 

Capturing the dynamics of the map can provide insights into 
the relationship between // and tj. To do this, we consider the fol - 
lowing symplectic map which approximates JUrmins kv 2008b), 
including the approximation ©, from one crossing of the SOS 
at time to to the next crossing at time t\ given by <I> : (to,Eo) — > 



£ l/2 
h/2 

t\ 
El 



. where 
= E 



a cos (to) +£>sin (to) 

( \ - 3 / 2 
to + a[-Ey 2 ) 

-3/2 



h/2- 



al-E 



1/2 



(9) 



E]j2 — a cos (t\) + bs'm (ty) 



and a, b and a are constants. The quantities t\ji and E\ji are the 
time and energy values of 7M3, respectively, when m-x, reaches a lo- 
cal maximum distance from the SOS with z(t\ii) = 0. The map is 
derived by approximating the change in energy of mj, on two time 
intervals in which we approximate its orbit by an orbit which es- 
capes parabolically. The first time interval (to, 1 1/2) corresponds to 
m3 moving away from the SOS. The second time interval (t\/2,t\) 
corresponds to the period in which mj returns to the SOS. It is 
clear that the change in energy is periodic in to and the trigono- 
metric terms in ([9} can be thought of as a lowest order Fourier ap- 
proximation to this change. The change in time is approximated 
by Keplerian motion over each time interval which means, for the 
chosen units, a = 7t/(2\/2). The map can be generalized as the it- 
erative map : (t„,E„) — » (t n +\,E n+ \) and the constants a and b 
are approximately proportional to e with 



a w 0.599 e/4 
b as 2.029 e/4. 

Sometimes it is more useful to write ((9) in the form 



t„ 
X„ 



= x, 



n-1 



■2a(-X n - 
4- 2£>sin(f„ 



-3/2 



(10) 



(11) 



for n = 1,2,3. 
generally X„ = 



-Ijl- 



3.1 Initial conditions 

Analogous to Moser's Do and D\ for {2j, we can define an open 
domain (Jo for which <!> is defined which is mapped into an open 
region U\ . Since time enters into the change in energy with period 
2jt, and we can transform energy values into velocity values by lO, 
we can consider Uq in polar coordinates where the angular argu- 
ment is determined by t and the radial argument is determined by 
v. An upper bound on allowable energy values in Uo is given by 

E i '(f) = -acos(f)-asin(f), (12) 

for t € [0,2n] which corresponds to E\ & = for which the map is 
undefined. All points in Uo get mapped to the open set U\ whose 




Figure 6. The boundaries BUq and dU\ for the regions £70 and Ui for e = 
0.61. 



boundary, dU\, is defined by 



E°(t) 



-acos (f) +&sin (t) 



(13) 



for t € [0,27t]. The boundaries dUo and 3(7 1 are depicted in Fig- 
ure [6] Allowable energy values in Uq at time t satisfy E < E? (t). 
To satisfy the physical constraints of the Sitnikov problem, energy 
values in the domain Uq are also bounded from below. From y), 
energy values on the SOS must satisfy 



1 

W) 



>0. 



(14) 



It has been shown bv lUrminskvl d2008bl) . that the dynamics of 
orbits with initial conditions in £7o for the map <5 are similar to the 
dynamics of orbits with initial conditions in Do for the map (]). More 
specifically, it was shown that th e map <£> satisfies lemmas similar 
to those proved bv lMoser! ( fT973l) which prove the existence of a set 
A £ Uo on which the dynamics are topologically equivalent to the 
shift map on the set of bi-infinite sequences. 

Initial conditions in Uo can be iterated forwards using © until 
the resulting orbits take on energy and time values which are out- 
side the domain Uq- A comparison of the Poinc are map and the 
approximate Poincare map <I> can be found in Urminskv (2008a) 
in which regions of initial values on the SOS are delineated by the 
number of excursions from the SOS an orbit makes before escap- 



where X = Ey 2 = E + acos(t ) + bsm(t ), and 3 _ 2 Definitions 



For a given orbit Z = {(tj,Ej)}f =0 computed by <S>, where N is the 
number of excursions from the SOS before 1113 escapes, the survival 
time is defined to be tj = f# — to- The growth of the logarithm of 
the solutions to the variational equations for the orbit Z is approxi- 
mated by, 



ln|8ZA. 



N 
i=0 



where Z^r = (t^,E^) and w,- is determined by, 
w,-_i 



Ji-i 



for i = l. 



.N. 



m-i 



(15) 



(16) 



in which |wrj[ = 1 is chosen at random and J,- is the lacobian of $ 
at time step i. After a few iterations w; is aligned with the unstable 
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Figure 7. The (ti,tj) scatter plot for 10000 uniformly distributed initial 
conditions in the domain of the map <I> for e = 0.61. The dashed line is 
the median curve associated with the scatter plot. 
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Figure 9. Median curves for the map <£> for varying eccentricity values. 
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Figure 8. Two part power law relationship for 100000 initial conditions for 
the map (e = 0.61) 
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Figure 10. The median curve on the interval containing 90 per cent of the 
points in Figure [8] On this interval a power law relationship with p ss 1 . 1 
fits the data best. 



direction associated with the solution at the rth time step. In a simi- 
lar way to equation {7}, we can express the relationship between 
and // as, 



1=0 

3.3 Numerical Results 

10000 random uniformly distributed initial conditions were cho- 
sen in Co such that, when iterated using equation {5), they escaped 
within 1000 iterations. If they failed to escape we did not include 
them in the calculations as there are regions of initial conditions in 
(Jo for which the corresponding orbits never escape. While calculat- 
ing the orbit we simultaneously compute < l 1 5b so as to determine the 
Lyapunov time by equation l |17t . Figure[7J shows the (f/,f<j) scatter 
plot where the dashed line represents the median curve associated 
with the scatter plot. Notice that the horizontal spread of the scatter 
plot for small time values found in Figure[3]is present in Figure[7j In 
addition, the density of the plotted points decreases as /; increases. 
In Figure[8]we plot the median curve for the (f; , tj) scatter plot 



for 100000 uniformly distributed initial conditions for e = 0.61. 
Again, there appears to be a two-part power law relationship for 
t\. The power law (TJ with p = 2.5 approximately fits the median 
curve on the interval 1 < ?/ < 9, whereas a power law with p = 1.07 
fits better for ?/ > 9. Figure [9] shows the median curves for various 
eccentricities of the binary. As in Figure[5] as the eccentricity of the 
binary increases, the two-part power law become apparent. 

Finally, we consider the time interval 4.4 <t\ < 129.1 which 
contains 90 per cent of the scatter points in Figure [7J This interval 
was chosen such that 5 per cent of the scatter points were in the 
region ?/ < 4.4 and 5 per cent of the points were in the region t\ > 
129.7. It was found that in this region a power law relationship with 
(3 ~ 1 . 1 best fitted the median curve (Figur efTOt. This contrasts with 
the result in Mikkola & Tanikawa] J2007h for the general 3-body 
problem where it was found that a power law with P ~ 1.8 roughly 
fits 90 per cent of the data on an interval .94 < ?/ < 35.2. Since 
the data in Figure [Tol is distributed over a larger interval, the data 
for larger t\ values dominates the approximation and the smaller 
power law approximation best fits the data. We shall demonstrate 
in section [4] that the distribution of for the map is different than 
that found for the general 3-body problem which may account for 
the discrepancy between the two results. 
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Figure 11. The boundary of the region B\. 



3.4 Delineating the region corresponding to escape after one 
excursion 

Initial conditions along the boundary dUo defined by U2\ lead to 
E\ji = 0. From equation (9), we can see that these initial conditions 
lead to a time fj/2 which is undefined. Initial conditions on E' (t) 
for t £ (0, Jt) are contained in the set D\ and so are in the domain of 
the inverse map 3> , Iterating these initial conditions backwards 
once gives the boundary of the region, B\, corresponding to penul- 
timate crossings of the SOS before escape. The boundary of B\ can 
be shown to be given parametrically by, 



U = r-2a(2i>sin(r))~ 3/2 
E- t = — 2fosin(f) — acos(ft) — fosin(f*). 



(18) 



This boundary is shown in Figure QT] using polar co-ordinates 
where the angular argument is time and the radial argument is the 
velocity of obtained from l[3j. Iterating initial values on this 
boundary forward in time using ((9} leads to energy values £3/2 = 
and from equation ^ we get, 

£3/2 =£i/2 + 2fesin(/i) =0. (19) 

Writing l |19l > in terms of £0 and to and rearranging gives, 

£() + a cos (to ) + b sin (to ) = 
-2isin^o + 2a(-£ , o-acos(/o)-isin(/ ))- 3 / 2 ) . 1201 

As the boundary of B\ spirals inside Uo (Figure [TH it approaches 
the boundary dUo- Consider Eg values on the boundary of B\ for a 
fixed to - Define 

e = 2a(-£ , o-acos(?o)-ksin(fo)r 3/ ' 2 (21) 

and note that — > °° as £0 — » E' (to). Using equation l !21b we can 
re-write J20b as 



-2/3 



= 2fcsm(f + e), (22) 

and as 6 — > °° solutions to J22b are asymptotically approximated by, 

Q + tO&kK, (23) 

for large iel Substituting equation i ]2 1 [ 1 into d23t we can obtain 
the approximation for the energy values on the boundary of B\ 



E^- 



( kit -tp 
\ 2a 



— a cos (to ) — b sin (to ) 



(24) 
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Figure 12. The (f;,^) scatter plot for orbits with initial values on the bound- 
ary of B\ for varying initial vector wq ( £ = 0.61). 



for large teZ. Initial conditions on B\ survive only one iteration 
of the map <&. Substituting l !24t into l[9} gives the survival time for 
orbits of initial conditions on the boundary of B\ as 

-3/2 



fd = fi-tow2<x 



(kit — to 
[ 2a 



-2/3 



-kK — tQ. 



(25) 



3.5 A functional relationship between t t and t d 

To derive a functional relationship between f/ and t d , we consider 
initial conditions on the boundary of B\ for a fixed to - For to = 7t/2, 
initial conditions on dB\ are approximated by, 



to 
Eq 



ji/2, 



(£- 1/2)71 \~ 2/3 



2a 



(26) 



for large teZ. Using the initial conditions i26\ and the correspond- 
ing survival time given by l !25t . we can compute the Lyapunov time 
ti from l !17t for a given initial |wo| = 1. For long lived orbits it may 
be assumed that wj will normally become aligned with the unsta- 
ble direction after only a few iterations. For initial conditions along 
3Bi on the other hand, this assumption is invalid since the orbit es- 
capes after only one excursion. The choice of wo has an important 
role in determining the growth of |5Zi| as it may not necessarily be 
aligned with the most unstable direction. 

The vector wo should be chosen so as to give the largest |5Zi| 
possible. Choosing wo in this way has the effect of making ?/ small. 
Figure [T2] shows the (f/,frf) scatter plot for initial values given by 
( 1 1 M for various choices of wo. It was found that the vector wo = 
[0, 1] maximizes In |8Zi | and ensures that In |8Zi | > for orbits of 
initial conditions on the boundary of B\ . 

For initial conditions in Si, the contribution of the denomina- 
tor in equation dl7t is ln[wi | since |wo| = 1. For initial conditions 
l !26t with initial vector wo = [0, 1] we can compute |wi | from equa- 
tion i Kit to obtain 



l + 3ta(-l) 



10/3 



(k- 



l/2)7t\ 5/3 



2a 



(27) 



For large k we have 
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Figure 13. Various approximations to the scatter (t/,t d ) plot of orbits which 
escape after one excursion from the SOS. 



Figure 14. The solid line is the median curve in Figure[8]for e = 0.61. The 
dashed line is the relationship 1351 with p = 0.16 and V = .12. 



and hence, for orbits of initial values J26t . 



(28) 



hi 



(34) 



ln|wi|~51n(fat)/3. (29) 

Substituting this result into dl7t . we find that the Lyapunov time for 
initial values on the boundary of B\ behaves like, 

3kn 



for large ieZ, 



(30) 



' 51n(ta) 

where the survival time is determined by J25l >. Since ~ kn we 
can rewrite OOt as 

3f rf 



(31) 



which for large t d is of order tj/\n(t,j). 

Figure [T3] shows various approximations to the scatter (tiJj) 
plot (represented by the solid line) for orbits which escape after one 
excursion from the SOS. A best fit for a curve in the form of the 
power law relationship |QJ is shown by the long dashed line where 
(3 = 1.07. This curve does not fit the scatter plot for small or large f/ 
very well. The dotted curve represents the function f/ = tj/ln | wj | , 
for |wj | given by \27\ , which approximates the scatter plot better 
than the power law relationship. Finally, the small dashed curve is 
the function // = 3/<//(51n(/j) which is only a slightly better ap- 
proximation to the scatter plot for small t d , but a much simpler one. 

3.6 A general functional relationship between ?/ and tj 

Now consider the orbits which escape after two excursions from the 
SOS. From equation (TTJ the Lyapunov time can be approximated 

by, 

| r^-n (32) 
ln|wi| +ln|w2| 

Following the result of the previous section, we assume that each 
term in the denominator is proportional to the log of the time spent 
between successive crossings of the SOS. For escape after two 
crossings these sum to tj and we take the two contributions to be 
ln(ft d ) and ln((l - f)t d ) where < / < 1, i.e. 

td 



Hftd)+H(i-f)td)' 

Re-arranging {33} we get 



(33) 



' Hf{i- f )tjy 

This relationship suggests fitting the data to curves which look like, 

td 



(35) 



ln(vf rf ) 

where p and v are constants. In fact, adding successive iterations 
produces a similar result so we take d35t as a general relationship 
and try to vary p and v to optimize a fit to the data. Figure[T4lshows 
the median curve (red solid curve) for 100000 initial values which 
escape within 1000 iterations of the map for e = 0.61. We were 
able to approximately fit the data with a function like i35l where 
p = 0.16andv = .12. 



4 DISTRIBUTIONS 

Now we look at the distributions of the various quantities studied 
so far and compare them to the results for the general 3-body prob- 
le m. In the numerical experime nts for the general 3-body problem 
in lMikkola & Tanikawal 1 120071) . the authors fit the marginal distri- 
bution of td (for large values of t d ) and large values of the ratio 
Z = tiijt\ to exponential probability density functions, finding that 



vfifd) ~ aexp(-a^), a= 1/250, 
V(Z) w pexp(-pZ), 0=1/45. 



(36) 
(37) 



The probability density for ?/ in the general 3-body problem was 
obtained numerically and a theoretical explanation was given in 
Mikk ola & Tanikawd ( 2007) as follows. The authors assume that 
t,j and Z are independent variables and so the marginal probability 
density function, \|/(f;), can be determined via 



v|/(f,) h j?>{ti-tdlZ)ysf{td)y{Z)dtddZ, 
= ap/(cw / + p) 2 . 



(38) 



but the 



Their numerical results do not quite match up with 
authors did fit their data with a function proportional to tf 2 . The 
problem may be that equation J3 8 b assumes that t d and Z are in- 
dependent variables. This does not seem to be likely considering 
Figure|3] It may be that, while the distributions provide satisfactory 
fits, an expression like l |38t may not always be true. 

Orbits computed with the map <5 also possess a similar power 
law relationship between and ?/, for large enough eccentricities, 
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DATA 
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ab/(at d +bf 
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Figure 15. The solid curve is the probability density of t d for the numerical 
results with £ = 0.61. The dashed curve was found to be the best approxi- 
mation of the stated form to the numerical results. 
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Figure 17. The solid curve is the probability density of Z for the numerical 
results with £ = 0.61. The dashed curve was found to be the best approxi- 
mation of the stated form to the numerical results. 
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Figure 16. The solid curve is the probability density of t d for the numerical 
results with £ = 0.61 for small t d . The dashed curve was found to be the 
best approximation of the stated form to the numerical results. 
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Figure 18. The solid curve shows the probability density of t\ for the nu- 
merical experiments with £ = 0.61. 



to those found in the general 3-body problem and the Sitnikov prob- 
lem. The map may provide some insights into the distributions for 
tj, Z and t[. Figure [J_5] shows the numerical results for the distri- 
bution of tj values for the map S> where e = 0.61. The t c [ values 
are binned into intervals of length 200 along the whole range of 
^ td ^ 70000 and the number of tj values in each bin is counted 
and divided by the total number of initial conditions. It was found 
that a probability density function 

>2 



y{t d )^a x bx/{ait d + biY 



(39) 



where a\ = 0.48 and b\ = 216.3 best represents the data in Fig- 
ure [T5] Figure [16] shows the results for < 1000 which demon- 
strates that this is a good approximation for small t d . Similarly, it 
was found that a probability density function for Z which fits the 
numerical data satisfactorily is given by 

^2 



x ¥{Z)^a 1 b 1 /{a 1 Z + b 1 y 



(40) 



where 02 = 10.35 and bi = 160.71 as shown in Figure [771 

The marginal probability density functions for tj and Z are 
quite different to those for for the general 3-body problem. This 
means that the predicted density function ( 138b may not necessarily 
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Figure 19. The solid curve shows the probability density of ?; for the nu- 
merical experiments with £ = 0.61 for small fy. 
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hold for the map. Proceeding as in d38t , assuming the variables t d 
and Z are independent, the density function for t[ can be determined 

by, 



«Pfo) = J 8(t, - t d /Z)V(t d )V{ tl )dt d dt, 



2b\a 2 \n f — — 
\a\ti 

{-a 2 b\ +b 2 a 1 ti) 3 

?b 2 

2a\b 2 ti\n — 

V"2, 



1+ln 



a\ti 



(-a 2 bi +b 2 a l ti) 2 
l+ln(^ 



(41) 



(— a 2 b\ + b 2 a\tiy {—a 2 b\+b 2 ci\ti) 2 

Figure Q~8] shows the numerical results for the map <!> for e = 0.61. 
The solid curve shows the data from the experiments, and the 
medium dashed curve is the predicted curve d4 1 b where g(f;) = 
Note that the curve does not quite fit the data from the exper- 
iments for small t[ (Figure[T9]l. The data was then fitted to a curve 
of the form, 



»P(f/) ~ ab/(at d +b) 2 , where a = 0.53, b = 2.75, 



(42) 



shown by the short dashed curve. Again, this failed to fit the data 
for small t\ . It was noted that the data is better approximated by the 



curve, 



f(ti)=32/tf. 



(43) 



Interestingly, this is a similar function to that which fits the numer- 
ical results for the general 3-body problem. This suggests that the 
probability distribution associated with the Lyapunov time may not 
necessarily be the result of the probability distribution of survival 
times and the probability distribution of the variable Z. 



t[ well. This may point to a more general property which may be 
valid for 3-body problems which experience large perturbations to 
an escaping mass. 
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5 CONCLUSIONS 

One of the results of this investigation is demonstrating a relation- 
ship between the Lyapunov time and the survival time for the Sit- 
nikov proble m which is similar to that for the general 3-body prob- 
lem found in iMikkola & Tanikawa (20071) . This is surprising as the 
Sitnikov problem is rather different compare d to the 3-body sys- 
tems discussed in Mi kkola & Tanikawal ( 120071) . With the use of an 
approximate Poincare map we were able to delineate regions of es- 
cape on a surface of section so as to construct initial conditions 
for the map. By studying the relationship between the Lyapunov 
time and the survival time with initial conditions in distinct escape 
regions, we were able to analytically obtain a new functional rela- 
tionship between t\ and t d given by t\ = pt d / ln(vt d ) where p and v 
are constants. As the (t[,t d ) scatter plots for the Sitnikov problem 
are similar to the (ti,t d ) scatter plots for the general 3-body prob- 
lem, we conjecture that the new functional relationship between t{ 
and t d presented above may also be valid for the general 3-body 
problem. 

Interestingly, the marginal distributions for the quantities t d 
and Z were found to be different than for the general 3-body prob- 
lem. The reason for this is not known although it may just be that in 
this study we considered the entirety of the numerical results and 
not just large values. Treating t d and Z as independent variables, 
we were able to derive a marginal distribution for ?/ which did not 
quite capture the numerical results. It seems unlikely though that t d 
and Z are independent quantities which may be why the theoretical 
predicted distribution *F(f/) did not quite match the numerical re- 
sults. Interestingly, as for the general 3-body problem, it was found 
that a function proportional to tf 2 fits the numerical distribution of 



